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Abstract 

Recently, a new family of integrators (Hamiltonian Boundary Value Meth- 
ods) has been introduced, which is able to precisely conserve the energy 
function of polynomial Hamiltonian systems and to provide a practical con- 
servation of the energy in the non-polynomial case. 

We settle the definition and the theory of such methods in a more general 
framework. Our aim is on the one hand to give account of their good behav- 
ior when applied to general Hamiltonian systems and, on the other hand, 
to find out what are the optimal formulae, in relation to the choice of the 
polynomial basis and of the distribution of the nodes. Such analysis is based 
upon the notion of extended collocation conditions and the definition of dis- 
crete line integral, and is carried out by looking at the limit of such family 
of methods as the number of the so called silent stages tends to infinity. 
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1 Introduction 



We consider canonical Hamiltonian problems in the form 

y = JVH{y), y(to) = yo e M^"^, (1) 

where J is a skew- symmetric constant matrix, and the Hamiltonian H{y) is as- 
sumed to be sufficiently differentiable. For its numerical integration, the problem 
is to find numerical methods which preserve H{y) along the discrete solution 
{yn}, since this property holds for the continuous solution y{t). 

So far, many attempts have been made inside the class of Runge-Kutta meth- 
ods, the most successful of them being that of imposing the symplecticity of the 
discrete map, considering that, for the continuous flow, symplecticity implies the 
conservation of H(y). Concerning symplectic integrators, a backward error anal- 
ysis permits to prove that they exactly conserve a modified Hamiltonian, even 
though this fact clearly does not always guarantee a proper qualitative behavior of 
the discrete orbits. 

On the other hand, it is possible to follow different approaches to derive geo- 
metric integrators which are energy-preserving. This has been done, for example, 
in the pioneering work [6], and later in [fT4ll . where discrete gradient methods are 
introduced and studied. An additional example of energy-preserving method is 
the Averaged Vector Field (AVF) method defined in [15J (see also [|4J). By the 
way, the latter method can be retrieved by the methods here studied. 

More recently, in O a new family of one-step methods has been introduced, 
capable of providing a numerical solution of ([T]), along which the energy 
function Hiy) is precisely conserved, in the case where this function is a polyno- 
mial (see also EolIIIlIIl). 

These methods, named Hamiltonian Boundary Value Methods {HBVMs here- 
after), may be also thought of as Runge-Kutta methods where the internal stages 
are split into two categories: 

- the fundamental stages, whose number, say s, is related to the order of the 
method; 

- the silent stages, whose number, say r, has to be suitably selected in order 
to assure the energy conservation property for a polynomial H(y) of given 
degree u; the higher is u, the higher must be r. 

The resulting method is denoted by HBVM(A;, s)Q where k = s + r is the total 
number of unknown stages. 

'The denomination HBVM with k steps and degree s was used in ||2l- 
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In [|2l [m it has also been shown that these new methods provide a practical 
conservation of the energy even in the non-polynomial case: the term "practical" 
means that, in many general situations, when the number of silent stages is high 
enough, the method makes no distinction between the function H{y) and its poly- 
nomial approximation, being the latter in a neighborhood of size e of the former, 
where e denotes the machine precision. 

Another relevant issue to be mentioned is that the computational cost for the 
solution of the associated nonlinear system is essentially independent of the num- 
ber of silent stages, and only depends on s (see [|2l[D)- This comes from the fact 
that the silent stages are actually linear combinations of the fundamental stages. 

These two aspects motivate the following question: what is, if any, the limit 
method when the number of silent stages grows to infinity? 



This question was first posed by Ernst Hairero who also provided a partial 
answer by stating formulae (|2TI) . which he called Energy Preserving variant of 
Collocation Methods {EPCMs, hereafter) [7J . We provide a proof of his statement 
by clarifying the connection between the limit formulae and HBVMs: we show 
that actually one can define several different limit methods^ each one associated 
to the specific polynomial basis, as well as to the choice of the abscissae distribu- 
tion, used to construct the sequence of HBVMs. For example, EPCMs are based 
upon the use of Lagrange polynomials, while, working with the shifted Legendre 
basis, yields to different limit methods, that we have called Infinity Hamiltonian 
Boundary Value Methods (in short, oo-HBVMs or HBVM(oo, s), being s the num- 
ber of the unknown fundamental stages). 

Our aim in this paper is threefold: 

1. We settle the definition of HBVMs in a more general framework, also de- 
riving the general formulation of the limit formulae 



In particular, we show that such limit coincides with EPCMs if the Lagrange 
polynomial basis is used (Section [2l). 

2. In Section[3l we introduce the new class of oo-HBVMs, which are the limit 
formulae corresponding to the HBVMs based upon the shifted Legendre 
polynomial basis. We prove that the order of such formulae is the same 
as the Gauss -Legendre methods, that is 2s (where s is the number of the 
unknown fundamental stages). 

^During the international conference "ICNAAM 2009", Rethymno, Crete, Greece, 18-22 
September 2009, after the talks, by the first two authors, where HBVMs were presented. 
^In the sense that they generate different discrete problems. 




lim HBVM(A;,s). 
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3. We mention the case where H{y) belongs to vector spaces different from 
that of polynomials, thus providing a natural (and trivial) generalization of 
the original formulae (see Section S]). Moreover, in the polynomial case, we 
determine the optimal distribution of the nodes (Section |5]). 

We stress that any finite approximation of EPCMs or oo-HBVMs based on 
quadratures leads back to HBVM(A;,s) methods, for k high enough. 

We address all the points listed above, by slightly modifying the approach 
followed to define the class of HBVMs in 01. 



2 Reformulation of Hamiltonian B VMs 

The key formula which HBVMs rely on, is the line integral and the related prop- 
erty of conservative vector fields: 

H{y^) - H{yo) = h [ &{to + ThfVH{a{U + r/i))dr, for any G M^-, 

(2) 

where a is any smooth function such that 



a 



(to) = 2/0, o-(io + h)= yi. (3) 



Here we consider the case where a{t) is a polynomial (of degree at most s), yield- 
ing an approximation to the true solution y{t) in the time interval [to, to + h]. The 
numerical approximation for the subsequent time-step, yi, is then defined by ([3]). 
After introducing a set of s distinct abscissae ci, . . . , Cs, (0 < q < 1)0 we set 

Yi = a(to + Cih), i = l,...,s, (4) 

so that cr(t) may be thought of as an interpolation polynomial, Yi, i = 1, . . . , s, 
being the internal stages. 

Let us consider the following expansions of cr(t) and a{t) for t G [to, to + Z^]: 

cr(to + Th) = 7i^j(^), o-(^o + Th) =yo + hy2 Ij / Pji^) da;, (5) 

where {-Pj(t)} is any suitable basis of the vector space of polynomials of degree 
at most s — 1 and the (vector) coefficients {7^} are to be determinedjl Before 
proceeding, one important remark is in order. 



'^As a convention, when c = is to be considered, as in the case of the Lobatto abscissae in 
[0, 1], then Co = is formally added to the abscissae ci, . . . , Cs, and the subsequent formulae are 
modified accordingly. 

^More general function spaces will be considered in the sequel. 
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Remark 1. As will be clear in a while, we observe that the numerical method 
which the following procedure will define is "basis-dependent" , in that to different 
choices of the basis {Pj{t)} there will, in general, correspond different numerical 
methods. In this section, in order to let the theory be presented as general as 
possible, we leave the basis not better specified. This will allow us to achieve the 
results listed at point 1. in the introduction. The question about how to choose 
the basis properly is faced in Section \3\ where oo-HBVMs will be introduced. 
Therefore, just in the present section, to avoid confusion, we will always specify 
what is the basis we are working with. This will be not necessary anymore starting 
from Section\3\ after determining the optimal basis. 

In this section we assume that H(y) is a polynomial, which implies that the 
integrand in Q is also a polynomial so that the line integral can be exactly com- 
puted by means of a suitable quadrature formula. It is easy to observe that in 
general, due to the high degree of the integrand function, such quadrature formula 
cannot be solely based upon the available abscissae {cj}: one needs to introduce 
an additional set of abscissae, ci, . . . , Cr, distinct from the nodes {q}, in order to 
make the quadrature formula exact: 

[ &{to + ThfVH{a{to + Th))dT= (6) 
Jo 

s r 

P,a{to + Cih)^VH{a{to + ah)) + J] ^i&{to + c,h)^VH{(x{to + ah)), 

i=l i=l 

where i = 1, . . . , s, and z = 1, . . . , r, denote the weights of the quadrature 
formula corresponding to the abscissae {q} and {q}, respectively, i.e.. 




(7) 




According to ifTTll . the right-hand side of Q is called discrete line integral, 
while the vectors 

Yi = a{tQ + Cih), i = l,...,r, (8) 

are called silent stages: they just serve to increase, as much as one likes, the 
degree of precision of the quadrature formula, but they are not to be regarded as 
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unknowns since, from ([5]), they can be expressed in terms of linear combinations 
of the fundamental stages (Hj). 

In [Q, the method HBVM(A;,s), with k = s + r is then defined by substituting 
the quantities in ^ into the right-hand side of ^ and by choosing the unknowns 
{7j} in order that the resulting expression vanishes. 

Instead of carrying out our computation on the right-hand side of Q, as was 
done in |[2|, we apply the procedure directly to the original line integral appearing 
in the left-hand side. Of course, since these two expressions are equal, the final 
formula will exactly match the HBVM(A;,s) method, written in a different guise. 

With this premise, by considering the first expansion in ([5]), the conservation 
property reads 

s „i 

J2lJ P,(r)ViJ((T(to + r/i))dr = 0, (9) 

which, as is easily checked, is certainly satisfied if we impose the following set of 
orthogonality conditions 

7i = Vj [ Pj{r)JVH{a{to + r/i))dr, j = l,...,s, (10) 

with {r]j} suitably nonzero scaling factors that will be defined in a while. Then, 
from the second relation of © we obtain, by introducing the operator 

LU-h)a{to + ch) = (11) 

a{to) + hy2vj Pj{x)dx / P,ir)fiaito + rh))dT, ce[0,l], 
j^i Jo Jo 

that a is the eigenfunction of L(JVH; h) relative to the eigenvalue A = 1: 

a = L{JVH;h)a. (12) 

Definition 1. Equation 4721) will be called the Master Functional Equation defin- 
ing o. 

Remark 2. 'We also observe that, from rfTOl) and the first relation in (|5l), one ob- 
tains the equations 



s „i 

d(to + c^h) = V VjPjici) / Pjir)JVH{a{to + Th))di 



(13) 

which may be viewed as extended collocation conditions according to / f77] Sec- 
tion 2], where the integrals are (exactly) replaced by discrete sums (see, e.g., ®- 

m). 
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To practically compute a, we set (see (H)) and dS])) 

s 

Yi = a{to + Cih) = yo + h^aij'jj, i = l,...,s, (14) 



where 



'0 



Pj{x)dx, i,j = l,...,s. (15) 



Inserting (flOl) into (fT4l) yields the final formulae which define the HBVMs class 
based upon the basis {Pj}: 

Y, = yQ + hj^ (^^r]ja^,P,{r)j JVH{a{to + rh))dr, i = 1, . . . , s. (16) 

The constants {r]j} have to be chosen in order to make the formula consistent. 
Problem ((T4l) - (fT6l) can be actually solved, provided that all the {r]j} are different 
from zerojj and the matrix 

/ /;'Pi(x)dx ... /;^Pi(x)dx \ 

\ f^^ Psix)dx ... Jq" Ps{x)dx J 

is nonsingular (which we shall obviously assume hereafter). Indeed, such a matrix 
allows us, by using ([5]), to reformulate equation (fT6l) in terms of the (unknown) 
fundamental stages dH). Let us now formally set 

f{y) = JVH{y), (17) 

and report a few examples for possible choices of the basis {Pj{x)}. 

1. In [lllil we have chosen {Pi(x), . . . , Ps{x)} as the Newton basis. This has 
allowed us the construction of a family of methods of order 2 and 4. 

2. In [|2l, the abscissae {cq = 0} U {q} U {q} are disposed according to a 
Lobatto distribution with k + 1 points in [0, ll and {Pi{x), . . . , Ps{x)} is 
the shifted Legendre basis in the interval [0, l]u Consequently, choosing in 
(fT6l) to = 0,h = 1, and /(y(r)) = Pjir), the consistency condition yields 

r]j = (^j^ P^{x)d:>^ = 2j - 1, J = 1, . . . , s, (18) 



^For example, the choice Pj{x) — x^^^, j = 1, . . . , s, would lead to r/i — 1, and rjj = 0, 
j = 2, . . . , s. This implies that, with this choice of the basis, a can only be a line (i.e., s = 1). 
'More precisely, Pj (x) is here the shifted Legendre polynomial of degree j — 1, j = 1, . . . , s. 
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which is exactly the value found in [|2l|. In such a case, it has been shown that 
the resulting method has order 2s, just the same as the generating Lobatto 
III A method (obtained for k = s). 

3. In a similar way, when using the Lagrange basis {ij{x)}, by setting /(i/(r)) = 
1, one obtains r]j = l/bj with 

bj = / ij{x)dx, £j{x) = Yl (19) 
Consequently, formulae (fT6l) become 

Y, = yo + hj^ ^^,{r)j JVH{a{to + r/i))dr, i = l,...,s. 

(20) 

Moreover, by introducing the new variables Ki = a{to + Cih), which are 
therefore related to the Yi as 

s 

Yi = yo + h^ttijKj, i = l,...,s, 
i=i 

system (|20|) can be recast in the equivalent form provided by the extended 
collocation conditions (fT3l) : 

1 

Ki = - l,{T)JWH{a{to + Th))dT, 2 = l,...,s. (21) 

Oi Jo 

Formulae (|2T]) are the ones E. Hairer proposed in the general case, that is for any 
kind of Hamiltonian function. They were called Energy Preserving variant of 
Collocation Methods (EPCMs) [7J. The above discussion then proves that if the 
integral can be substituted by a finite sum, as in the case where H{y) is a poly- 
nomial, formulae ([161) . and consequently (|2TI) . become a HBVM(A;, s), with a 
suitable value of A;|j 

For sake of completeness, we report the nonlinear system associated with the 
HBVM(A;, s) method, in terms of the fundamental stages {Yj} and the silent stages 
{Yj} (see dS])), by using the notation (fTTl) . In this context, they represent the dis- 
crete counterpart of (fT6l) . and may be directly retrieved by evaluating, for example, 

**The same clearly happens when the integral is only approximated by a finite sum. 
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the integrals in (fT6l) by means of the (exact) quadrature formula introduced in 

s j s \ r I s \ 

Yi = yo + h Y^^nYl VjaijPjici) f{Yi) + ^ A VjaijPj{ci) f{Yi) 
1=1 \j=i J 1=1 \j=i J 

(22) 

s / s r \ 

= yo + hJ2vjaiAY,PiPjici)f{Yi) + Y,PiPj{ci)f{Yi)\, i = l,...,s. 
j=i \i=i 1=1 ) 

From the above discussion it is clear that, in the non-polynomial case, supposing 
to choose the abscissae {q} so that the sums in (l22l) converge to an integral as 
r = A; — s — )■ oo, the resulting formula is (fT6l)FI Consequently, EPCMs may be 
viewed as the limit of HBVMs family, when the Lagrange basis is considered, as 
the number of silent stages grows to infinity. 

The above arguments also imply that HBVMs may be as well applied in the 
non-polynomial case since, in finite precision arithmetic, HBVMs are indistin- 
guishable from their limit formulae (fT6l) . when a sufficient number of silent stages 
is introduced. The aspect of having a practical exact integral, for k large enough, 
was already stressed in [|2l [T0l[Tni . 



3 Infinity Hamiltonian Boundary Value Methods 

As is easily argued (and emphasized in Remark [B, the choice of the basis along 
which cr(to + tK) is expanded (see ([5])), somehow influences the shape of the final 
formulae (fT6l) . that is, to two different polynomial bases there may correspond 
two different families of formulae^ The question then naturally arises about the 
best possible choice of the basis to consider. This issue has been a crucial point in 
devising the class of HBVMs in [2] and deserves a particular attention^ 

Indeed, we recall that our final goal is to devise methods that make the sum ([91), 
representing the line integral, vanish. This is accomplished by the orthogonality 
conditions (flOl) . whose effect is to make null each term of the sum in It follows 
that such conditions are in general too demanding, in that they are sufficient but 
not necessary to get conservativeness. In fact, the sum could in principle vanish 
even in the case when two ore more of its terms are different from zero. This extra 
constraint may affect the general properties of the conservative methods we are 
interested in, and in particular their order. 

This was a problem already encountered in IfTTll where the authors realized 
that the use of the Newton basis didn't assure the expected growth of the order of 

^This obvious requirement for the abscissae will be always assumed in the sequel. 
'°For example, see the method presented in subsection l4.2l 

"The argument presented here is the analog of the one appearing in yjj Remark 3.1]. 
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the resulting method when the degree of the polynomial a (to + Th) was increased. 
This barrier has been definitively overcome in [2], where it was understood that 
the proper polynomial basis to be used by default was that of the shifted Legendre 
polynomials in the interval [0, 1]. We emphasize that, contrary to what happens 
for the Lagrange and Newton bases, the Legendre polynomials are orthogonal and 
symmetric in the interval [0, 1] and in addition they are abscissae-free, that is they 
by no means depend on the specific distribution of the abscissae {cj} adopted. 
This, in turn, implies that the Master Functional Equation (fT2l) is independent of 
the choice of both the abscissae {cj} and {q}: the only requirement being that ^ 
holds truelll 

From the above arguments, it is clear that the orthogonality conditions (flOl) . 
i.e., the fulfillment of the Master Functional Equation (fT2l) . is only a sufficient 
condition for the conservation property ^ to hold. Such a condition becomes 
also necessary, when the basis {Pj} is orthogonal. 

Theorem 1. Let {Pj} be an orthogonal basis on the interval [0, 1]. Then, assum- 
ing H{y) to be suitably differentiable, ([P]) implies that each term in the sum has to 
vanish. 

Proof Let us consider, for simplicity, the case of an orthonormal basis, and 
the expansion 

g{r) = VH{a{to + rh)) = ^ p,P,{t), p, = (P,, g), £ > 1, 
where, in general, 

U\9)= I f{r)9{r)dT. 
Jo 

Substituting into yields 

i=i j=i V i>i / j=i 

Since this has to hold whatever the choice of the function H(y), one concludes 
that 

7jPi = 0, j = l,...,s.D (23) 

Remark 3. In the case where {Pj} is the shifted Legendre basis, from f [2?]) one 
derives that 

jj = Spj, i = l,...,s, 

where S is any nonsingular skew -symmetric matrix. The natural choice S = J 
then leads to 4701) . with rjj = {Pj, Pj)^^ = 2j — 1. 

^^We emphasize that this is not the case when using, for example, the Lagrange basis. 
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The use of the Legendre basis allows the resulting methods to have the best 
order and stability properties that one can expect. This aspect is elucidated in 
the two theorems and the corollary below, which represent the main result of the 
present work. 

Although, up to now, we have maintained the treatment of HBVMs at a general 
level, it is clear that, in view of the result presented in Theorem [2l when the curve 
cr(to + Th) is assumed of polynomial type, we will implicitly adopt the Legendre 
basisl^ This important assumption will be incorporated in the HBVM methods 
from now on: if needed, the use of any other kind of basis will be explicitly stated, 
in order not to create confusion. 

Taking into account the consistency conditions (fTSl) . formula (fT6l)-(fT5l) takes 
the form: 

Yi=yo + hj^ (^^{2j -l)ai,P,{T)j JVH{(r{to + rh))dT, z = l,...,s. 

(24) 

If the Hamiltonian H{y) is a polynomial, the integral appearing at the right- 
hand side is exactly computed by a quadrature formula, thus resulting into a 
HBVM(A;,s) method with a sufficient number of silent stages. As already stressed 
in the previous section, in the non-polynomial case such formulae represent the 
limit of the sequence HBVM(A;,s), as /c — )■ oo. 

Definition 2. We call the new limit formula dH]) an Infinity Hamiltonian Boundary 
Value Method (in short, oo-HBVM or HBVM{oo, s) ). 

We emphasize that, in the non-polynomial case, (l24l) becomes an operative 
method, only after that a suitable strategy to approximate the integral is taken 
into account (see the next section for additional examples). In the present case, 
if one discretizes the Master Functional Equation (fTTI) - (fT2l) . HBVM(A;, s) are 
then obtained, essentially by extending the discrete problem (l22l) also to the silent 
stages dS]). In order to simplify the exposition, we shall use (fTTl) and introduce the 
following notation: 

{ti} = {ci}U{ci}, {w,} = {A}U{A}, yi = a{to+tih), h = f{(T{to+tih)). 

(25) 

The discrete problem defining the HBVM(/c, s) then becomes, with r]j = 2j — 1, 



yi = yo + h^r]j / Pj{x)dx ^ ujePj{ti)fi, i = l,...,k. (26) 
j=i -^0 e=i 



'^Actually, the term Hamiltonian Boundary Value Method has been coined in f2\, after intro- 
ducing the Legendre basis. 



11 



We can cast the set of equations in vector form, by introducing the vectors y = 
{yf, vlY, e = (1, . . . , 1)^ G Ml", and the matrices X, P G Ml'''', with 

lij = r Pj{x)Ax, Vij=P,{ti), (27) 
A = diag(?7i,...,?7s), = diag(a;i, . . . , w^), 

as 

y = e®yo + hiXAV^n) ® / /(y), (28) 

with an obvious meaning of f{y). Consequently, the method can be seen as a 
Runge-Kutta method with the following Butcher tableau: 



tl 








tk 






Ui ... Uk 



(29) 



Remark 4. We observe that, provided that the matrix A is independent of the basic 
abscissae {q} (as in the case of the Legendre basis), the role of such abscissae 
and of the silent abscissae {q} is interchangeable. This is not true, for example, 
for the Newton and Lagrange bases. 

The following result then holds true. 

Theorem 2. Provided that the quadrature has order at least 2s (i.e., it is exact for 
polynomials of degree at least 2s — Ij, HBVM(k,s) has order p = 2s = 2 deg(cr), 
whatever the choice of the abscissae Ci, . . . , c^. 

Proof From the classical result of Butcher (see, e.g., |l9l Theorem 7.4]), the 
thesis follows if the simplifying assumptions C(s), B{p), p > 2s, and D{s — 1) 
are satisfied. By looking at the method (|28l) -(l29l). one has that the first two (i.e., 
C{s) and B{p), p > 2s) are obviously fulfilled: the former by the definition of 
the method, the second by hypothesis. The proof is then completed, if we prove 
D{s — 1). Such condition can be cast in matrix form, by introducing the vector 
e = (1, . . . , 1)"^ G W^^, and the matrices 

g = diag(l,...,s-l), D = diag(ti,...,tfc), = (tr^) G M^^^-\ 

(see also (ITT]) ) as 

QV^n {lAV'^n) = (ee^ - V^D) n, 

i.e.. 
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Vhl^nVQ = (e - DV) . (30) 
Since the quadrature is exact for polynomial of degree 2s — 1. one has 




Pi(x)xMx^ , z = l,...,s, j = l,...,s-l, 

where the last equality is obtained by integrating by parts, with Sn the Kronecker 
symbol. Consequently, 



i = l,...,k, j = l,...,s-l, 
that is, (l30l) . where the last equality follows from the fact that 

s „i 

J2vePeit) / Peix)xMx = t\ 



j = i,...,s-i.n 



Concerning the stability, the following result holds true. 

Theorem 3. For all k such that the quadrature formula has order at least 2s = 
2deg(cT), HBVM(k,s) is perfectly A-stable, whatever the choice of the abscissae 
Ci , . . . , Cg. 

Proof As it has been previously observed, a HBVM(A;, s) is fully charac- 
terized by the corresponding polynomial a which, for k sufficiently large (i.e., 
assuming that ^ holds true), satisfies the Master Functional Equation (fTTI) - (fT2l) . 
which is independent of the choice of the nodes Ci, . . . , Cg (since we consider 
the Legendre basis). When, in place of f{y) = J'VH{y) we put the test equa- 
tion f(y) = \y, we have that the collocation polynomial of the Gauss-Legendre 
method of order 2s, say as, satisfies the Master Functional Equation, since the in- 
tegrands appearing in it are polynomials of degree at most 2s — 1, so that cr = cr^. 
The proof completes by considering that Gauss-Legendre methods are perfectly 
A-stable. □ 

A worthwhile consequence of Theorems [2] and [3] is that one can transfer to 
HBVM(cxD, s) all those properties of HBVM(A;,s) which are satisfied starting from 
a given k > ko on: for example, the order and stability properties. 
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Corollary 1. Whatever the choice of the abscissae Ci, . . . , Cg, HBVM{oo, s) (|24l) 
has order 2s and is perfectly A-stable. 

Remark 5. From the result of CorollaryU} it follows that HBVM{oo, s) has order 
2s and is perfectly A-stable for any choice for the abscissae ci, . . . , Cs. Since such 
abscissae can be arbitrarily chosen, we can formally place them at the roots of 
the Gauss -Legendre polynomial of degree s. On the other hand, by considering 
that, at such abscissae, by setting {ii{c)} and {bi} the corresponding Lagrange 
polynomials and quadrature weights, respectively ( see f [79l) ), 



1 1 ^ 

— / Pj{x)ii{x)dx = —y2brPj{Cr)iiiCr) = Pj{Ci), 
Jo 

one obtains ( with rjj = 2j — 1 and by using the notation rfTTI) ): 

s „i 

a'{to + Cih) = Y,TljPj{ci) / P,(r)/(a(to + r/i))dr 



T)f{a{to + Th))dr 



II' (j2vjPAr) jyAxMx)dx^ f{a{to + Th))dT 



h 



1 



- / ii{T)f{cr{to + Th))dT, z = 



bi JO 



Consequently, for any choice of the abscissae {q}, HBVM{oo, sYprovide the same 
polynomial a as the "optimal EPCMs" ([27]) of order 2s /77o Conversely, an 
EPCM is optimal (i.e., it has order 2s) only when the abscissae Ci, . . . , define 
a quadrature formula of order at least 2s — 1, whereas different choices result in 
methods of lower order Theorem 1 ]. 

Remark 6. We also observe that, due to the choice of the shifted Legendre poly- 
nomial basis (see d24l) ) 



HBVM{oo, s) = Um HBVM{k, s 



A;— >oo 



whatever is the choice of the fundamental abscissae {q}. Consequently, for all k 
large enough, so that the Master Functional Equation 4721) holds true (e.g., in the 
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In this sense, they are equivalent, even though they generate different discrete problems. 
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case of a polynomial Hamiltonian H{y)), all HBVM{k, s) provide the same poly- 
nomial a of degree s, independently of the choice of the abscissae {q}. Hence, 
they are equivalent to each other. This result doesn 't change in the case where 
H{y) is not a polynomial, provided that H{y) is sufficiently differentiable. In this 
case, in fact, one formally obtains, in place of the Master Functional Equation 
(fT2l) . an equation of the form 

ak = L{JVH;h)ak + Mh), (31) 

where ipk{h) = 0{h'^''~^^'^), qu being the degree of precision of the quadrature at 
the right-hand side in ©, so that ^ oo as k oo. From ([72]) and ( 1371) . one 
then obtains that as h ^ 0, assuming that f is Lipschitzian with constant fi, and 
for a suitable constant M independent of h: 

\Wk - <y\\ < hfiM\\ak - cr\\ + ||^fc(/i)||, 

i.e., 

\Wk - < (1 - hnM)-^\\Mh)\\ = 0{h'''-'+^) ^0, k^oo. 

One then concludes that, when using finite precision arithmetic, ak is indistin- 
guishable from a, for all k large enough. 

Example 1. As previously mentioned, for the methods studied in [2], based on 
a Lobatto distribution of the nodes {cq = 0, ci, . . . , Cs} U {ci, . . . , Ck-s}, one 
has that deg((T) = s, so that the order of HBVM(k,s) turns out to be 2s, with a 
quadrature satisfying B{2k). 

Example 2. For the same reason, when one considers a Gauss distribution for 
the abscissae {ci, . . . , Cg] U {ci, . . . , Ck-s], one also obtains a method of order 2s 
with a quadrature satisfying B{2k). This case will be further studied in Section\5\ 

Remark 7. Finally, we also mention that, from Remark^ HBVM(k,s) are sym- 
metric methods]^ provided that the abscissae {ti} (see ( l25l) ) are symmetrically 
distributed (see also /|2j/j. 

4 Generalization of Hamiltonian B VMs 

The approach that has allowed the construction of methods that conserve energy 
functions of polynomial type is quite general: that is, by no means it depends 
on the particular vector space generating the curve a{t) nor on the quadrature 

According to the time reversal symmetry condition defined in [S] p. 218]. 
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technique used. As was emphasized in IfTTl Section 2], it solely relies on the 
following two ingredients: the definition of discrete line integral and the extended 
collocation conditions (fT3l) . which zero the line integral 

Therefore, in a more general context, this procedure can be formalized as fol- 
lows. One first picks a curve cr(to + T/i), T G [0, 1], joining two points of the phase 
space = (j{tQ) and yi = (T(to + h)- Such a curve is assumed to lie in a proper 
finite dimensional vector space W = span{Pi(x), . . . , Ps{x)}, where now Pj{x), 
j = 1, . . . , s, are any linearly independent functions. Therefore the curves a(t) 
and a(t) will admit an expansion in the form ([5]). 

The fundamental hypothesis, for this approach to work, is that the choice of 
W must guarantee that the functions Pj(T)V H(a{tQ + rh)) appearing in ^ (and 
hence ^{tYV H{a{t))) be elementary integrable, that is they are required to admit 
a primitive that can be expressed in terms of elementary functions. If this is the 
case, all the steps performed to obtain (fT6l) may be repeated with the integral 
substituted by the primitive. 

This represents a generalization of what done for polynomial Hamiltonian 
functions not only because the vector space W may be generated by non-poly- 
nomial functions but also because the analytic solution of the line integral may be 
carried out by any available technique. Hereafter, we report a couple of examples 
in the class (I24l)n 



4.1 A method of order two 

We consider a separable Hamiltonian function (for simplicity we assume m = 1) 



H{q,p) = V{p)-U{q). 
Let a{t) be the segment joining = (go,Po)^ to yi = 

cr(to + rh) = yo + r{yi - yo). 
We have Cq = 0, Ci = 1, and the corresponding method (l24l) becomes 

\ 



/ 9i - go \ 




h 




Pi -Po 




\ h } 





\ 



y\p^ + ^(Pi -po))cir 
+ ^(gi - go))dr 







P\ -Po 

f/(gi) - 



qi - go 



(32) 



/ V{pi) - V{p,) \ 



(33) 



'^While the method in Section|4T|is equivalently obtainable by applying either ( l24l i or (fTSl l. the 
same is not true for the fourth-order method derived in Section l4~2l where. the use of the Lagrange 
basis, would produce a coefficient 62 = (appearing as a denominator in the resulting formulae 
(ED). 
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Formula (l33l) is one of the simplest discrete gradient methods due to Itoh and 
Abe lfT2l . whose general form, for non-separable Hamiltonian functions with one 
degree of freedom, reads 



/ qi 


- qo 








h 




= J 


Pi 


-Po 






\ 


h 


I 





qi - qo 

H{qi,Pi) - H{qi,po) 



\ 



(34) 



Pi -Po 



J 



The vector appearing at the right-hand side of (|34l) is obtained by replacing the 
partial derivatives of H(q, p) with increments along the q and p axes. Method (l34l) 
is in general first order and nonsymmetric. However, when confined to separable 
Hamiltonian systems, it turns out to be second order and symmetric0 



4.2 A method of order four 



To construct a method of order four in the form (|24l) applied to (l32l) . we pick a 
curve a{t) of degree two, based upon the abscissae cq = 0, ci = 1/2, and C2 = 1. 
Such a method has been already described in [fTTll for polynomial Hamiltonian 
functions: here we consider its generalization to the non-polynomial case. Setting 
Yi = {qi/2,Pi/2)^ and, observing that Y2 = iqi,Pi)^, the two components of the 
curve a {to + rh) are 



(Tiito + rh) 
(T2(to + rh) 



2{qo - 2gi/2 + qi)r'^ - (3go - 4gi/2 + qi)r + go 
2{po - 2pi/2 + Pi)r2 - (3po - 4pi/2 + Pi)r + po 



Consequently, ((24l) becomes 



Yi 



(35) 





+ h 



j\-\r + \)V\a,{to + rh))dr 



\-\r + -^U\a,{to + rh))<\r 



\ 



Yo 





\ 



V\a2{to + rh))Ar 



U'{ai{to + rh))Ar 



\ 



J 



(36) 

Substituting (l35l) into (l36l) we obtain a system in the unknowns qi/2, P1/2, qi. Pi- 
Looking at (l36l) . we realize that even in the simpler case of a system deriving 



'^A generalization of ( |34] | introduced in [13] also becomes method ( |33] | when applied to Hamil- 
tonian functions in the form ( |32] |. 
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from a Hamiltonian function in the form (l32l) . the elementary integrability of the 
integrals in (l24l) is not a priori guaranteed. This means that, in this case, we cannot 
arrive at a general formula analogous to (l33l) . in terms of U{q) and V{p). 

On the other hand, in several cases of interest, such primitive can be explicitly 
computed: hereafter we report a significant example, which we shall use later in 
the numerical tests in Section [5j 

Example 3. The role of this example is also to show that, when finite precision 
arithmetic is used, it may be not convenient to use the infinite version of the meth- 
ods, even if the integrals can be analytically evaluated. This will be evident from 
the numerical results in Section l54l The system we consider is the one defined by 
the Hamiltonian function 

H{q,p) = a{\ogq-q) + h{\ogp-p), (37) 

where a and b are positive constants. The associated system ([T]) reads 

q= h(--l\, p = -a(--l\. (38) 



kP J \<i 
This system is strictly related to the Lotka-Volterra model 

q= bq{l-p), 
p = -ap{l - g), 



(39) 



in that system (|39l) may be recast as the Poisson system y = ^J^^ JVH{y), where 

ri{q,p) = — ^ is called integrating factor. 

Systems |39I) and (1381) share the same Hamiltonian function (|37]) as first in- 
tegral and, consequently, they share the same curves as trajectories in the phase 
plane. Method (l36l) applied to (l39l) reads 



/ gl/2 - go \ 

h/2 

Pl/2 -PO 

\ h/2 J 



-b + 
■ ^arctanh(' 



3l log(|po/pi|) 



1 b Po- 



3l/2 + 7pi 



-3po 



4pi/2-Pi 



1 ' 2 Ci PO 

) — arctanh( 



2pi/2+Pl 
Po-4pi/2+3pi 



y • (arctanh(' 



3^Jog(_|go/£i|)_ 
4 9o-2(ji/2+gi 



-3qo+4qi/2-qi 

C2 



1 a go-8gi/2 + 7gi 

2 C2 <3o-2i3i/2+<3i 
^ 90 -4£i/2+3gi 



(40) 



) — arctanh(- 



C2 



/ 


91 


-go \ 




( 








h 1 










Pi 


-Pa 








V 




h ) 




-+'cT2 



/ 6 _ 2fc rarctanh( ^"'"P^f+'^^ ) + arctanh(^^^^i^%i^)j \ 



arctanh(*^i%^) + arctanh( '^°''^y^+^^ ) 



(41) 
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where 

Ci = (pI + lQpli2 + pI- 8poPi/2 - 2poPi - 8pl/2Pl)^/^ 
C2 = (go + 16g^^2 + - 2gogi - ^qi/2qif'^- 

5 HBVMs based upon Gauss quadrature 

As anticipated in Example [21 we now study the properties of the HBVM(/c, s) 
which is defined over the set of k distinct abscissae, 

{ti, . . . , 4} = {ci, . . . , cj U {ci, . . . , Cfc_4, 

coinciding with the Gauss-Legendre nodes in [0, 1], i.e., the roots of the shifted 
Legendre polynomial of degree k. The corresponding polynomial a has then de- 
gree s. By virtue of Theorems |2] and [3] (see also Remark |7]), such methods are 
symmetric, perfectly A-stable, and of order 2s. They reduce to Gauss-Legendre 
collocation methods, when k = s, and are exact for polynomial Hamiltonian func- 
tions of degree provided that 

k>Y- (42) 

By recalling what stated in Remark [6l for all k sufficiently large so that ^ holds, 
HBVM(A;, s) based on the k Gauss-Legendre abscissae in [0, 1] are equivalent 
to HBVM(A;, s) based on A; + 1 Lobatto abscissae in [0, 1] (see [|2l), since both 
methods define the same polynomial a of degree 50 

As matter of fact, we have run HBVM(A;, s) based on Gauss-Legendre nodes, 
and HBVM(A;, s) based on the Lobatto nodes, obtaining the same results on the 
polynomial test problems reported in |l2l, which are briefly recalled in the sequel. 

5.1 Test problem 1 

Let us consider the problem characterized by the polynomial Hamiltonian (4.1) in 
0, 

Hip,g) = P-^^- + ^ + '--'-+^-, (43) 
^■^'^^ 3 2 30 4 3 6' ^ ^ 

having degree z/ = 6, starting at the initial point yo = (q'(0),p(0))^ = (0, 1)^, 
so that H(yo) = 0. For such a problem, in (5] it has been experienced a numer- 
ical drift in the discrete Hamiltonian, when using the fourth-order Lobatto IIIA 
method with stepsize h = 0.16, as confirmed by the plot in Figure [T] When us- 
ing the fourth-order Gauss-Legendre method the drift disappears, even though the 



**In the non-polynomial case, they converge to the same HBVM(oo, s), as /c — ?> cx). 
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Table 1 : Maximum difference between the numerical solutions obtained through 
the fourth-order HBVM(A;, 2) methods based on Lobatto abscissae and Gauss- 
Legendre abscissae for increasing values of k, problem (l45l) . 10^ steps with step- 
size h = 0.1. 



k 


h = 0.1 


2 


3.97-10-1 


4 


2.29 • 10"^ 


6 


2.01 ■ 10-*^ 


8 


1.37- 10-" 


10 


5.88 • 10-13 



Hamiltonian is not exactly preserved along the discrete solution, as is shown by 
the plot in Figure [21 On the other hand, by using the fourth-order HBVM(6,2) 
with the same stepsize, the Hamiltonian turns out to be preserved up to machine 
precision, as shown in Figure [3l since such method exactly preserves polynomial 
Hamiltonians of degree up to 6. In such a case, the numerical solutions obtained 
by using the Lobatto nodes {cq = 0, ci, . . . , cg = 1} or the Gauss-Legendre nodes 
{ci, . . . , cg} are the same. 

5.2 Test problem 2 

The second test problem, having a highly oscillating solution, is the Fermi-Pasta- 
Ulam problem (see |[8l Section 1.5.1]), defined by the Hamiltonian 

H{p,q) = ^^{pli-i+pli) + Y^(^2i- q2i~if + ^iq2i+i- q2i)^ , (44) 

i=l i=l i=0 

with go = q2m+i = 0, m = 3, = 50, and starting vector 

Pi = 0, qi = {i-l)/10, i = l,...,6. 

In such a case, the Hamiltonian function is a polynomial of degree 4, so that 
the fourth-order HBVM(4,2) method, either when using the Lobatto nodes or the 
Gauss-Legendre nodes, is able to exactly preserve the Hamiltonian, as confirmed 
by the plot in Figure[6l obtained with stepsize h = 0.05. Conversely, by using the 
same stepsize, both the fourth-order Lobatto IIIA and Gauss-Legendre methods 
provide only an approximate conservation of the Hamiltonian, as shown in the 
plots in Figures m and [5l respectively. 
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Figure 1: Fourth-order Lobatto IIIA method, h = 0.16, problem (l43l) : drift in the 
Hamiltonian. 
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Figure 2: Fourth-order Gauss-Legendre method, h = 0.16, problem (|43l) : H 
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Figure 3: Fourth-order HBVM(6,2) method, h = 0.16, problem dH: H ^ 10"^^ 
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Figure 4: Fourth-order Lobatto IIIA method, h = 0.05, problem dH): \H-Ho\ 
10-3. 
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Figure 5: Fourth-order Gauss-Legendre method, h 
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Figure 6: Fourth-order HBVM(4,2) method, h = 0.05, problem dH): \H - Hq] 

10-1^ 



5.3 Test problem 3 (non-polynomial Hamiltonian) 

In the previous examples, the Hamikonian function was a polynomial. Neverthe- 
less, as observed above, also in this case HBVM(A;,s) are expected to produce a 
practical conservation of the energy when applied to systems defined by a non- 
polynomial Hamiltonian function that can be locally well approximated by a poly- 
nomial. As an example, we consider the motion of a charged particle in a magnetic 
field with Biot-Savart potentiallll It is defined by the Hamiltonian [2J 

(45) 



with g = a/x^ + y"^, a = e Bq, m is the particle mass, e is its charge, and Bq is 
the magnetic field intensity. We have used the values 

m = 1, e = —1, Bq = 1, 

with starting point 

x = 0.5, y = lO, z = 0, x = -0.1, y = -0.3, i = 0. 

By using the fourth-order Lobatto IIIA method, with stepsize h = 0.1, a drift is 
again experienced in the numerical solution, as is shown in Figure |71 By using 
the fourth-order Gauss-Legendre method with the same stepsize, the drift disap- 
pears even though, as shown in Figure [H the value of the Hamiltonian is pre- 
served within an error of the order of 10^'^. On the other hand, when using the 
HBVM(6,2) method with the same stepsize, the error in the Hamiltonian decreases 
to an order of 10^^^ (see Figure |9l), thus giving a practical conservation. Finally, 
in Table [T] we list the maximum absolute difference between the numerical so- 
lutions over 10^ integration steps, computed by the HBVM(/c, 2) methods based 
on Lobatto abscissae and on Gauss-Legendre abscissae, as k grows, with stepsize 
h = 0.1. As expected, the difference tends to 0, as k increases, since the two 
sequences of methods tend to the same limit, given by the HBVM(oo, 2) (see ((24l) 
with s = 2). 

5.4 Test problem 4 (non-polynomial Hamiltonian) 

We finally solve the Hamiltonian system (l38l) by using the Itho-Abe method (1331) . 
the fourth-order formula (l40l)-(l4T]). and the HBVM(10,2), which has order four 

This kind of motion causes the well known phenomenon of aurora borealis. 
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Figure 7: Fourth-order Lobatto IIIA method, h = 0.1, problem (|45l) : drift in the 
Hamiltonian. 
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Figure 8: Fourth-order Gauss-Legendre method, h 
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Figure 9: Fourth-order HBVM(6,2) method, h = 0.1, problem (|45]): \H - Ho\ 

10-15. 



and degree of precision 10 (that is, according to (l42l) . it precisely conserves the 
energy of polynomial Hamiltonians of degree up to 10). We have set a = 6 = 1 in 
formula (|37|) . and integrated over a time interval [0, 5000] with stepsize h = 0.5 
and {qo,Po) = (0.5,0.5) as initial condition. 

Figure[lO]reports the numerical Hamiltonian function associated with the three 
methods. The occurrence of jumps in the first two graphs (left picture) is due to the 
fact that both formulae (l33l) and (|401)-(I4T]) may become ill-conditioned for certain 
values of the state vector. For example (see Figure [TT]), at the two consecutive 
times t = 2830.5 and t = 2831, the state vectors associated with the Itoh-Abe 
method (|33l) are, respectively, 

[qi,pi] ~ (0.39988668, 1.4216560)^, [^2,^2] ^ (0.39988872, 0.67130503)^, 

which shows that qi may be very close to q2 even for large values of h. This 
causes some cancellation in the subtraction at the right-hand side of (l33l) and, 
hence, a jump of the subsequent branch of the numerical trajectory on a different 
level curve. However, since, in general, the numerical trajectory densely fills the 
level curve H{y) = H^yo), it may be argued that the occurrence of such jumps 
are systematic and frequent when the dynamics is traced over a long time. The 
use of finite arithmetic eventually destroys the theoretical conservation property. 
A similar argument may be applied to discuss the behavior of the fourth-order 
method (|40l)-(|4T]). 

Although the HBVM method does not provide a theoretical conservation of 
the energy, as is the case for the above cited methods, its behavior in finite arith- 
metic would suggest the opposite (see the right picture in Figure [TOl) . as already 
emphasized at the beginning of Example [3l 

6 Conclusions 

In this paper, the newly introduced Hamiltonian Boundary Value Methods (HB- 
VMs), a class of numerical methods able to exactly preserve polynomial Hamilto- 
nians of any degree, have been re-derived in a unifying framework. Such frame- 
work relies on the use of line integrals, which are approximated by suitable dis- 
crete counterparts (actually, they are exact, when the Hamiltonian is a polyno- 
mial). In this context, the limit of the methods, as the number of the so called 
silent stages tends to infinity, is easily obtained. When the underlying polyno- 
mial basis upon which the HBVM is constructed is the Lagrange basis, such limit 
formulae coincide with the recently introduced Energy Preserving variant of Col- 
location Methods; if instead one uses the shifted Legendre polynomial basis, the 
corresponding HBVMs have the highest possible order and so do their limit for- 
mulae, the Infinity Hamiltonian Boundary Value Methods, independently of the 
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Figure 10: Left picture: absolute error of the Hamiltonian function (ITT]) evalu- 
ated along the numerical solutions computed by the Itoh-Abe method (l33l) (lower 
curve) and formula (l40l)- (|4T]) (upper curve). The jumps are symptomatic of ill- 
conditioning of the formulae for certain values of the solution. Right picture: the 
same kind of plot produced by the HBVM formula of order 4 and A; = 10 Gaussian 
abscissae (\H - Ho\ ^ lO^^^^^ 
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Figure 11: Trajectory in the phase plane computed by the Itoh-Abe method (1331) . 
The small circles locate the solution at the two consecutive times t = 2830.5 and 
t = 2831. The very close values of the variable q for such two points causes loss 
of significant digits in the subsequent branch of the trajectory. 
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considered abscissae. Any limit formula, when discretized, fall into the HBVMs 
class. Possible extensions of the approach have been also sketched, along with a 
number of numerical tests. Such tests confirm that, in the limit of the silent stages 
tending to infinity, all HBVMs with s (unknown) fundamental stages tend to the 
same limit method, which is characterized by the eigenfunction (relative to the 
unit eigenvalue) of a certain operator, which is independent of the choice of the 
abscissae. 
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